#include <iostream>
#include <stdexcept>
#include <vector>
#include <ginac/ginac.h>
#include "sector_decomposition/sector_decomposition.h"

int main()
{

  try{

    using namespace sector_decomposition;
    using namespace GiNaC;

    std::cout << "testsuite sector decomposition" << std::endl;
    std::cout << std::endl;

    // --------------------------------------------------------------
    int verbose_level = 1;

    CHOICE_STRATEGY = STRATEGY_C;

    // --------------------------------------------------------------
    monte_carlo_parameters mc_parameters = monte_carlo_parameters( 5, 15, 10000, 100000 );

    // --------------------------------------------------------------

    symbol eps("eps");
    int order = -4;

    int n     = 7;
    int loops = 2;

    symbol minus_sinv("minus_sinv"), minus_tinv("minus_tinv"), minus_uinv("minus_uinv");
    symbol minus_p12("minus_p12"), minus_p22("minus_p22"), minus_p32("minus_p32"), minus_p42("minus_p42");
    symbol m12("m12"), m22("m22"), m32("m32"), m42("m42"), m52("m52"), m62("m62"), m72("m72");

    std::vector<ex> parameters;
    symbol x1("x1"), x2("x2"), x3("x3"), x4("x4"), x5("x5"), x6("x6"), x7("x7");
    parameters.push_back(x1);
    parameters.push_back(x2);
    parameters.push_back(x3);
    parameters.push_back(x4);
    parameters.push_back(x5);
    parameters.push_back(x6);
    parameters.push_back(x7);

    integration_data global_data = integration_data(parameters, eps, order);

    std::vector<exponent> nu_minus_1(n,exponent(0,0));

    // eq.5 of hep-ph/0305234
    ex U = (x1+x2+x3)*(x4+x5+x6+x7) + (x4+x5)*(x6+x7);
    ex F = (x2*x3*(x4+x5+x6+x7)+x2*x4*x6+x3*x5*x7)*minus_sinv // typo x3 x5 x6 --> x3 x5 x7
      +x1*x4*x7*minus_tinv
      +x1*x5*x6*minus_uinv
      +x1*(x5*x7+x2*(x4+x5+x6+x7))*minus_p12
      +x1*(x4*x6+x3*(x4+x5+x6+x7))*minus_p22
      +(x6*x7*(x1+x2+x3+x4+x5)+x3*x4*x7+x2*x5*x6)*minus_p32
      +(x4*x5*(x1+x2+x3+x6+x7)+x3*x5*x6+x2*x4*x7)*minus_p42
      + x1*m12 + x2*m22 + x3*m32 + x4*m42 + x5*m52 + x6*m62 + x7*m72;

    F = F.subs( minus_sinv == 1 );
    F = F.subs( minus_tinv == 1 );
    F = F.subs( minus_uinv == 1 );
    F = F.subs( minus_p12 == 0 );
    F = F.subs( minus_p22 == 0 );
    F = F.subs( minus_p32 == 0 );
    F = F.subs( minus_p42 == 0 );
    F = F.subs( m12 == 0 );
    F = F.subs( m22 == 0 );
    F = F.subs( m32 == 0 );
    F = F.subs( m42 == 0 );
    F = F.subs( m52 == 0 );
    F = F.subs( m62 == 0 );
    F = F.subs( m72 == 0 );

    std::vector<ex> poly_list;
    poly_list.push_back(U);
    poly_list.push_back(F);

    std::vector<exponent> c(poly_list.size());
    c[0] = exponent( n-(loops+1)*2, loops+1 );
    c[1] = exponent( -n+2*loops, -loops );
    for (int k=0; k<n; k++) 
      {
	c[0].sum_up(nu_minus_1[k]);
	c[1].subtract_off(nu_minus_1[k]);
      }

    integrand my_integrand = integrand(nu_minus_1, poly_list, c);

    // --------------------------------------------------------------

    monte_carlo_result res = do_sector_decomposition(global_data, my_integrand, mc_parameters, verbose_level);

    // --------------------------------------------------------------

  } catch (std::exception &e)
    {
      std::cout << "Exception : " << e.what() << std::endl;
    }

  return 0;
}

